1 Week 05 : 23 and 25 September 2008

We now begin the discussion of software issues such as step size control (Sects. continuous output ([HNW93] Sect. II.5) and automatic stiffness detection ([HW96] Sect. IV.2) schemes for embedded Runge–Kutta methods. 2 Stepsize control Suppose we take a code such as rksuite and try to reverse engineer it to extract the RAS of the methods it implements. Consider running it on the equation y ′ = λy and seeing when u n →0; the code requires an error tolerance ǫ, not too large, and varies h to provide a solution with local error ǫ in each step. When hλ is outside the RAS, the solution will be inaccurate and the code will reduce h to maintain accuracy. Thus the numerical solution will always →0 for ℜλ < 0, and the determination of the RAS is a little more subtle. The RAS could be determined by finding the cost of the method for each λ: in the RAS, the step size is O(ǫ −1/p) as the error ǫ→0, while outside the RAS the step size is determined by stiffness to be O(1/|λ|). This is tricky to decide, so we need to know a little more about how such codes choose the step size to control the error. The goal of step size selection strategies is to produce a solution which is very likely to be within a user–specified tolerance ǫ of the exact solution, at minimal cost (necessarily depending on ǫ like O(ǫ −1/p) for an order–p method). The classical step–doubling approach based on Richardson extrap-olation is instructive and will be useful, but modern codes use an alternate approach called " embedded " Runge–Kutta methods. It is important to control the error because engineering computations are worth a lot more when they include error estimates or bounds. The overall approach to controlling the error of a one–step method of order p is straightforward. Suppose we just took a step of size H from u n−1 to u n and computed an error estimate E n , and now we want to determine the largest step size h that will give us u n+1 with error estimate E n+1 ≤ ǫ. The step h may also be required to satisfy other constraints such as: a singularity occurs at internal time t * , where we would like to stop and restart, or the user has requested output at certain times T j. …